# Replication Archive for: 
# Coppock, Alexander and Donald P. Green. 2020. 
# "Do Belief Systems Exhibit Dynamic Constraint?" 
# The Journal of Politics, Forthcoming.

rm(list = ls())
library(tidyverse)
library(estimatr)
library(rsample)
library(rmeta)

study_1_mturk <- read_rds("data/clean/study_1_mturk_cleaned.rds")
study_1_elite <- read_rds("data/clean/study_1_elite_cleaned.rds")
study_1_lucid <- read_rds("data/clean/studies_1_2_replications_lucid_cleaned.rds")

study_1_mturk <- study_1_mturk %>% filter(Z != "climate")
study_1_lucid <- study_1_lucid %>% filter(Z != "climate")

dv_df <- read_rds("data/raw/dv_df.rds")
dv_df_mturk <- filter(dv_df, dv_cat != "climate")
dv_df_elite <- filter(dv_df, dv_cat != "climate")
dv_df_lucid <- filter(dv_df, dv_cat != "climate")

# Bootstrapping -----------------------------------------------------------


pwa <- function(data, dvs){
  
  all_fits <-
    dvs %>%
    map( ~ formula(paste0(., "_w1 ~ Z"))) %>%
    map( ~ lm_robust(., data = data)) %>%
    map(tidy) %>%
    bind_rows 
  
  all_fits <-
    all_fits %>%
    mutate(dvs = gsub(outcome, pattern = "_w1", replacement = "")) %>%
    filter(term != "(Intercept)") %>%
    left_join(dv_df) %>%
    mutate(treatment = factor(term,
                              levels = c("Zpaul", "Zwallstreet", "Zveterans", "Zamtrak"),
                              labels = c("Treatment: Flat Tax", "Treatment: Wall Street", "Treatment: Veterans", "Treatment: Amtrak")),
           dv_cat_label = factor(dv_cat,
                                 levels = c("flat", "wall", "vets","amtrak"),
                                 labels = c("Outcomes: Flat Tax", "Outcomes: Wall Street", "Outcomes: Veterans", "Outcomes: Amtrak")),
           dv_labels_rev = factor(dv_labels, levels = rev(levels(dv_labels))),
           dv_cat_harmony = str_split_fixed(pattern = ": ", as.character(dv_cat_label), n = 2)[,2],
           treatment_harmony = str_split_fixed(pattern = ": ", as.character(treatment), n = 2)[,2],
           same_cat = if_else(dv_cat_harmony == treatment_harmony, "Within Domain", "Across Domain"))
  
  ests <-
    all_fits %>%
    group_by(same_cat) %>%
    summarize(pwa = weighted.mean(estimate, w = (1/std.error^2))) %>%
    spread(same_cat, pwa) %>%
    mutate(ratio = `Across Domain` / `Within Domain`)
  
  return(ests)
}


# sims <- 1000
sims <- 10

mturk_resamples <- bootstraps(study_1_mturk, times = sims)
elite_resamples <- bootstraps(study_1_elite, times = sims)
lucid_resamples <- bootstraps(study_1_lucid, times = sims)

suppressMessages({
  suppressWarnings({
    mturk_bootstraps <-
      map_df(mturk_resamples$splits, ~ (pwa(., dvs = dv_df_mturk$dvs)))
    elite_bootstraps <-
      map_df(elite_resamples$splits, ~ (pwa(., dvs = dv_df_mturk$dvs)))
    lucid_bootstraps <-
      map_df(lucid_resamples$splits, ~ (pwa(., dvs = dv_df_mturk$dvs)))
  })
})


mturk_est <- study_1_mturk %>% pwa(dvs = dv_df_mturk$dvs)
mturk_uncertainty <- 
mturk_bootstraps %>%
  summarize_all(.funs = list(
    sd = ~ sd(.),
    quantile025 = ~ quantile(., probs = c(0.025)),
    quantile975 = ~ quantile(., probs = c(0.9755))
  )) %>%
  pivot_longer(cols = everything())
  

elite_est <- study_1_elite %>% pwa(dvs = dv_df_elite$dvs)
elite_uncertainty <- 
  elite_bootstraps %>%
  summarize_all(.funs = list(
    sd = ~ sd(.),
    quantile025 = ~ quantile(., probs = c(0.025)),
    quantile975 = ~ quantile(., probs = c(0.9755))
  )) %>%
  pivot_longer(cols = everything())

lucid_est <- study_1_lucid %>% pwa(dvs = dv_df_lucid$dvs)
lucid_uncertainty <- 
  lucid_bootstraps %>%
  summarize_all(.funs = list(
    sd = ~ sd(.),
    quantile025 = ~ quantile(., probs = c(0.025)),
    quantile975 = ~ quantile(., probs = c(0.9755))
  )) %>%
  pivot_longer(cols = everything())





meta.summaries(
  d = as.numeric(c(mturk_est[1, 3], elite_est[1, 3], lucid_est[1, 3])), 
  se = as.numeric(c(mturk_uncertainty[3,2], elite_uncertainty[3,2], lucid_uncertainty[3,2]))
)


mturk_est
mturk_uncertainty

elite_est
elite_uncertainty

lucid_est
lucid_uncertainty
